Análisis de Expresión Diferencial I

Licenciatura en Ciencias Genómicas, Transcriptómica, Semestre 2025-2

1. Motivación

El análisis de expresión diferencial identifica genes o moléculas con cambios significativos en sus niveles de expresión bajo diferentes condiciones, como enfermedades o tratamientos. Es crucial para comprender mecanismos biológicos, descubrir biomarcadores, desarrollar terapias dirigidas y avanzar en la medicina personalizada. En bioinformática, el DEA permite procesar y analizar grandes volúmenes de datos de tecnologías como RNA-Seq o microarrays, ayudando a identificar patrones clave en redes biológicas. Su importancia radica en generar conocimiento sobre enfermedades, guiar investigaciones y acelerar el desarrollo de herramientas diagnósticas y terapéuticas.

1.1 Setup de datos

Vamos a usar este set de datos: https://www.ncbi.nlm.nih.gov/geo/download/?acc=GSE63310&format=file

## Download data
url <- "https://www.ncbi.nlm.nih.gov/geo/download/?acc=GSE63310&format=file"
utils::download.file(url = url,
                     destfile = "GSE63310.tar",
                     mode = "wb")
utils::untar(tarfile = "GSE63310.tar",
             exdir = file.path(tempdir(),"GSE63310"))

## Unzip files
files <- c("GSM1545535_10_6_5_11.txt", "GSM1545536_9_6_5_11.txt", "GSM1545538_purep53.txt",
  "GSM1545539_JMS8-2.txt", "GSM1545540_JMS8-3.txt", "GSM1545541_JMS8-4.txt",
  "GSM1545542_JMS8-5.txt", "GSM1545544_JMS9-P7c.txt", "GSM1545545_JMS9-P8c.txt")
for(i in paste(files, ".gz", sep=""))
  R.utils::gunzip(file.path(tempdir(),"GSE63310/",i), overwrite=TRUE)

2. EdgeR / Limma

EdgeR es una librería de R que se usa para realizar análisis de expresión diferencial de genes provenientes de datos de RNA-seq. Esta librería emplea métodos estadísticos para experimentos multigrupo basados en modelos lineares generalizados (generalized linear models: glms) los cuales son eficientes para experimentos multifactoriales independientemente de su complejidad. También implementa métodos Bayesianos que permiten la estimación de la variación biológica genética-específica, incluso para experimentos con los mínimos niveles de replicación biológica.

• Documentación (EdgeR)
• DOI: https://doi.org/10.1093/bioinformatics/btp616

if (!require("BiocManager", quietly = TRUE))
    install.packages("BiocManager", quiet = TRUE)

if (!require(edgeR, quietly = TRUE))
    BiocManager::install("edgeR",quiet = TRUE)

Limma es una librería de R que fue creada para el ánalisis de expresión diferencial para datos provenientes de microarreglos. Usa principalmente modelos lineales para experimentos multigrupos.

• Documentación (limma)
• DOI: https://doi.org/10.1093/nar/gkv007

if (!require("BiocManager", quietly = TRUE))
    install.packages("BiocManager", quiet = TRUE)

if (!require(limma, quietly = TRUE))
    BiocManager::install("limma",quiet = TRUE)

2.1 Load libraries

suppressPackageStartupMessages(library(limma))
suppressPackageStartupMessages(library(Glimma))
suppressPackageStartupMessages(library(edgeR))
suppressPackageStartupMessages(library(Mus.musculus))

2.2 Read files

x <- readDGE(file.path(tempdir(),"GSE63310/",files),
             columns = c(1,3))
class(x)

## [1] "DGEList"
## attr(,"package")
## [1] "edgeR"

dim(x)

## [1] 27179     9

2.2.1 Reorganize information

## Assign sample type and sample names
samplenames <- substring(colnames(x),67)
samplenames

## [1] "10_6_5_11" "9_6_5_11"  "purep53"   "JMS8-2"    "JMS8-3"    "JMS8-4"   
## [7] "JMS8-5"    "JMS9-P7c"  "JMS9-P8c"

colnames(x) <- samplenames
group <- as.factor(c("LP", "ML", "Basal", "Basal", "ML", "LP",
                     "Basal", "ML", "LP"))
x$samples$group <- group
lane <- as.factor(rep(c("L004","L006","L008"), c(3,4,2)))
x$samples$lane <- lane
x$samples

2.2.2 Add gene names

genes <- AnnotationDbi::select(Mus.musculus,
                        keys = rownames(x),
                        columns = c("SYMBOL","TXCHROM"),
                        keytype = "ENTREZID")

## Remove duplicates
genes <- genes[!duplicated(genes$ENTREZID),]
head(genes)

## Add genes
x$genes <- genes
x

## An object of class "DGEList"
## $samples
##                                                                                           files
## 10_6_5_11 C:\\Users\\chomb\\AppData\\Local\\Temp\\RtmpwZFj5a/GSE63310//GSM1545535_10_6_5_11.txt
## 9_6_5_11   C:\\Users\\chomb\\AppData\\Local\\Temp\\RtmpwZFj5a/GSE63310//GSM1545536_9_6_5_11.txt
## purep53     C:\\Users\\chomb\\AppData\\Local\\Temp\\RtmpwZFj5a/GSE63310//GSM1545538_purep53.txt
## JMS8-2       C:\\Users\\chomb\\AppData\\Local\\Temp\\RtmpwZFj5a/GSE63310//GSM1545539_JMS8-2.txt
## JMS8-3       C:\\Users\\chomb\\AppData\\Local\\Temp\\RtmpwZFj5a/GSE63310//GSM1545540_JMS8-3.txt
## JMS8-4       C:\\Users\\chomb\\AppData\\Local\\Temp\\RtmpwZFj5a/GSE63310//GSM1545541_JMS8-4.txt
## JMS8-5       C:\\Users\\chomb\\AppData\\Local\\Temp\\RtmpwZFj5a/GSE63310//GSM1545542_JMS8-5.txt
## JMS9-P7c   C:\\Users\\chomb\\AppData\\Local\\Temp\\RtmpwZFj5a/GSE63310//GSM1545544_JMS9-P7c.txt
## JMS9-P8c   C:\\Users\\chomb\\AppData\\Local\\Temp\\RtmpwZFj5a/GSE63310//GSM1545545_JMS9-P8c.txt
##           group lib.size norm.factors lane
## 10_6_5_11    LP 32863052            1 L004
## 9_6_5_11     ML 35335491            1 L004
## purep53   Basal 57160817            1 L004
## JMS8-2    Basal 51368625            1 L006
## JMS8-3       ML 75795034            1 L006
## JMS8-4       LP 60517657            1 L006
## JMS8-5    Basal 55086324            1 L006
## JMS9-P7c     ML 21311068            1 L008
## JMS9-P8c     LP 19958838            1 L008
## 
## $counts
##            Samples
## Tags        10_6_5_11 9_6_5_11 purep53 JMS8-2 JMS8-3 JMS8-4 JMS8-5 JMS9-P7c
##   497097            1        2     342    526      3      3    535        2
##   100503874         0        0       5      6      0      0      5        0
##   100038431         0        0       0      0      0      0      1        0
##   19888             0        1       0      0     17      2      0        1
##   20671             1        1      76     40     33     14     98       18
##            Samples
## Tags        JMS9-P8c
##   497097           0
##   100503874        0
##   100038431        0
##   19888            0
##   20671            8
## 27174 more rows ...
## 
## $genes
##    ENTREZID  SYMBOL TXCHROM
## 1    497097    Xkr4    chr1
## 2 100503874 Gm19938    <NA>
## 3 100038431 Gm10568    <NA>
## 4     19888     Rp1    chr1
## 5     20671   Sox17    chr1
## 27174 more rows ...

2.3 Data preprocessing

## counts per million
cpm.obj <- cpm(x)

## log counts per million transformation
lcpm.obj <- cpm(x, log = TRUE)

## Mean and median
L <- mean(x$samples$lib.size) * 1e-6
M <- median(x$samples$lib.size) * 1e-6
c(L, M)

## [1] 45.48855 51.36862

## The average library size is about 45.5 million and the minimum
## log-CPM should be log2(2/45.5) = -4.51 which means that
## a count of 0 maps to a log-CPM value of -4.51
summary(lcpm.obj)

##    10_6_5_11          9_6_5_11          purep53             JMS8-2       
##  Min.   :-4.5074   Min.   :-4.5074   Min.   :-4.50743   Min.   :-4.5074  
##  1st Qu.:-4.5074   1st Qu.:-4.5074   1st Qu.:-4.50743   1st Qu.:-4.5074  
##  Median :-0.6847   Median :-0.3589   Median :-0.09513   Median :-0.0901  
##  Mean   : 0.1714   Mean   : 0.3312   Mean   : 0.43559   Mean   : 0.4089  
##  3rd Qu.: 4.2913   3rd Qu.: 4.5601   3rd Qu.: 4.60081   3rd Qu.: 4.5475  
##  Max.   :14.7632   Max.   :13.4952   Max.   :12.95700   Max.   :12.8513  
##      JMS8-3            JMS8-4            JMS8-5            JMS9-P7c      
##  Min.   :-4.5074   Min.   :-4.5074   Min.   :-4.50743   Min.   :-4.5074  
##  1st Qu.:-4.5074   1st Qu.:-4.5074   1st Qu.:-4.50743   1st Qu.:-4.5074  
##  Median :-0.4281   Median :-0.4064   Median :-0.07152   Median :-0.1704  
##  Mean   : 0.3225   Mean   : 0.2529   Mean   : 0.40428   Mean   : 0.3708  
##  3rd Qu.: 4.5772   3rd Qu.: 4.3199   3rd Qu.: 4.42513   3rd Qu.: 4.6031  
##  Max.   :12.9578   Max.   :14.8520   Max.   :13.19491   Max.   :12.9413  
##     JMS9-P8c      
##  Min.   :-4.5074  
##  1st Qu.:-4.5074  
##  Median :-0.3300  
##  Mean   : 0.2749  
##  3rd Qu.: 4.4355  
##  Max.   :14.0102

2.4 Remove genes that are lowly expressed

## Remove all the genes that have 0 counts across the nine samples
table(rowSums(x$counts==0)==9)

## 
## FALSE  TRUE 
## 22026  5153

## Filter by expression
keep.exprs <- filterByExpr(x,group = group)
x <- x[keep.exprs,,keep.lib.sizes=FALSE]
dim(x)

## [1] 16624     9

2.5 Normalize gene expression

Se necesita normalizar debido a que los datos de RNA-Seq contienen sesgos técnicos y biológicos que pueden distorsionar los resultados si no se corrigen. Las razones principales para normalizar son: diferencias de profundidad de secuenciación, diferencias en la composición génica, comparabilidad entre muestras, impacto de la longitud de los genes, preparación para métodos estadísticos.

Métodos comunes de normalización:

• TMM (Trimmed Mean of M-values): Corrige la profundidad de secuenciación y la composición génica.
• RPKM/FPKM (Reads/Fragments Per Kilobase Million): Normaliza por la longitud del gen y el tamaño total de las muestras.
• TPM (Transcripts Per Million): Similar a RPKM, pero mejora la comparabilidad entre muestras.
• DESeq2/EdgeR Scaling Factors: Ajustan los conteos para análisis estadísticos robustos.

## Normalize
x <- calcNormFactors(x,method = "TMM") # Trimmed mean of M-values
x$samples$norm.factors

## [1] 0.8943956 1.0250186 1.0459005 1.0458455 1.0162707 0.9217132 0.9961959
## [8] 1.0861026 0.9839203

Comparación entre normalizar o no

x2 <- x
x2$samples$norm.factors <- 1
x2$counts[,1] <- ceiling(x2$counts[,1]*0.05)
x2$counts[,2] <- x2$counts[,2]*5

Visualización:

col <- c("#A6CEE3","#1F78B4","#B2DF8A","#33A02C","#FB9A99","#E31A1C","#FDBF6F","#FF7F00","#CAB2D6")
par(mfrow=c(1,2))
lcpm <- cpm(x2, log=TRUE)
boxplot(lcpm, las=2, col=col, main="")
title(main="A. Example: Unnormalised data",ylab="Log-cpm")
x2 <- calcNormFactors(x2)  
x2$samples$norm.factors

## [1] 0.05770899 6.08287835 1.22023972 1.16478991 1.19661094 1.04659233 1.15048074
## [8] 1.25431164 1.10901983

lcpm <- cpm(x2, log=TRUE)
boxplot(lcpm, las=2, col=col, main="")
title(main="B. Example: Normalised data",ylab="Log-cpm")

2.6 Unsupervised clustering of samples

• Usar una MultiDimensional Scaling plot (MDS): esta muestra la similitudes y disimilitudes entre las muestras con algoritmo no supervisado. La primera dimensión representa el fold-change que mejor separa a las muestras y que explica la mayor proporción de la variancia de los datos. Aquí nos podemos dar cuenta si se deben de hacer algunas correcciones del batch.

lcpm <- cpm(x, log=TRUE)

## Define colors
par(mfrow=c(1,2))
col.group <- group
levels(col.group) <-  RColorBrewer::brewer.pal(nlevels(col.group), "Set1")
col.group <- as.character(col.group)

## Define colors
col.lane <- lane
levels(col.lane) <-  RColorBrewer::brewer.pal(nlevels(col.lane), "Set2")
col.lane <- as.character(col.lane)

## Plot MDS
plotMDS(lcpm, labels=group, col=col.group)
title(main="A. Sample groups")
plotMDS(lcpm, labels=lane, col=col.lane, dim=c(3,4))
title(main="B. Sequencing lanes")

2.7 Differential expression analysis

2.7.1 Create a matrix

Genes que están expresados en diferentes niveles. Los modelos lineales se ajustan a los datos bajo la asunción de que los datos siguen una distribución normal (lo cual no es cierto porque siguen una distribución binomial negativa o una de poisson).

Este diseño experimental quita el intercepto del primer factor group pero un intercepto permanece con el segundo factor lane. Es esencial entender cómo interpretar los coeficientes en el modelo.

design <- model.matrix(~0 + group + lane)
colnames(design) <- gsub("group","",colnames(design))
design

##   Basal LP ML laneL006 laneL008
## 1     0  1  0        0        0
## 2     0  0  1        0        0
## 3     1  0  0        0        0
## 4     1  0  0        1        0
## 5     0  0  1        1        0
## 6     0  1  0        1        0
## 7     1  0  0        1        0
## 8     0  0  1        0        1
## 9     0  1  0        0        1
## attr(,"assign")
## [1] 1 1 1 2 2
## attr(,"contrasts")
## attr(,"contrasts")$group
## [1] "contr.treatment"
## 
## attr(,"contrasts")$lane
## [1] "contr.treatment"

contr.matrix <- makeContrasts(
   BasalvsLP = Basal-LP, 
   BasalvsML = Basal - ML, 
   LPvsML = LP - ML, 
   levels = colnames(design))
contr.matrix

##           Contrasts
## Levels     BasalvsLP BasalvsML LPvsML
##   Basal            1         1      0
##   LP              -1         0      1
##   ML               0        -1     -1
##   laneL006         0         0      0
##   laneL008         0         0      0

2.7.2 Removing heteroscedascity from count data

Para los datos de RNA-seq la varianza no es independiente de la media sin importar si los datos sufren una transformación logarítmica o no. Los modelos que asumen una distribución binomial negativa asumen una relación cuadrática media-varianza.

El plot de voom ilustra la relación media-varianza. Este plot muestra una trend que decrece entre la media y la varianza como resultado de la combinación de la varianza técnica en los experimentos de secuenciación y de la varianza biológica entre las réplicas de las muestras de las diferentes poblaciones celulares.

Experimentos con una variación biológica alta resultan en tendencias más planas.

Experimentos con una variación biológica baja reusltan en tendencias decrecientes pronunciadas.

Si no se tienen más de 3 réplicas se opta por usar un GLM o si se desea usar un modelo específico para RNA-seq.

par(mfrow=c(1,2))
v <- voom(x, design, plot=TRUE)
v

## An object of class "EList"
## $genes
##   ENTREZID SYMBOL TXCHROM
## 1   497097   Xkr4    chr1
## 5    20671  Sox17    chr1
## 6    27395 Mrpl15    chr1
## 7    18777 Lypla1    chr1
## 9    21399  Tcea1    chr1
## 16619 more rows ...
## 
## $targets
##                                                                                           files
## 10_6_5_11 C:\\Users\\chomb\\AppData\\Local\\Temp\\RtmpwZFj5a/GSE63310//GSM1545535_10_6_5_11.txt
## 9_6_5_11   C:\\Users\\chomb\\AppData\\Local\\Temp\\RtmpwZFj5a/GSE63310//GSM1545536_9_6_5_11.txt
## purep53     C:\\Users\\chomb\\AppData\\Local\\Temp\\RtmpwZFj5a/GSE63310//GSM1545538_purep53.txt
## JMS8-2       C:\\Users\\chomb\\AppData\\Local\\Temp\\RtmpwZFj5a/GSE63310//GSM1545539_JMS8-2.txt
## JMS8-3       C:\\Users\\chomb\\AppData\\Local\\Temp\\RtmpwZFj5a/GSE63310//GSM1545540_JMS8-3.txt
## JMS8-4       C:\\Users\\chomb\\AppData\\Local\\Temp\\RtmpwZFj5a/GSE63310//GSM1545541_JMS8-4.txt
## JMS8-5       C:\\Users\\chomb\\AppData\\Local\\Temp\\RtmpwZFj5a/GSE63310//GSM1545542_JMS8-5.txt
## JMS9-P7c   C:\\Users\\chomb\\AppData\\Local\\Temp\\RtmpwZFj5a/GSE63310//GSM1545544_JMS9-P7c.txt
## JMS9-P8c   C:\\Users\\chomb\\AppData\\Local\\Temp\\RtmpwZFj5a/GSE63310//GSM1545545_JMS9-P8c.txt
##           group lib.size norm.factors lane
## 10_6_5_11    LP 29387429    0.8943956 L004
## 9_6_5_11     ML 36212498    1.0250186 L004
## purep53   Basal 59771061    1.0459005 L004
## JMS8-2    Basal 53711278    1.0458455 L006
## JMS8-3       ML 77015912    1.0162707 L006
## JMS8-4       LP 55769890    0.9217132 L006
## JMS8-5    Basal 54863512    0.9961959 L006
## JMS9-P7c     ML 23139691    1.0861026 L008
## JMS9-P8c     LP 19634459    0.9839203 L008
## 
## $E
##         Samples
## Tags     10_6_5_11  9_6_5_11   purep53     JMS8-2    JMS8-3    JMS8-4    JMS8-5
##   497097 -4.292165 -3.856488 2.5185849  3.2931366 -4.459730 -3.994060 3.2869677
##   20671  -4.292165 -4.593453 0.3560126 -0.4073032 -1.200995 -1.943434 0.8442767
##   27395   3.876089  4.413107 4.5170045  4.5617546  4.344401  3.786363 3.8990635
##   18777   4.708774  5.571872 5.3964008  5.1623650  5.649355  5.081611 5.0602470
##   21399   4.785541  4.754537 5.3703795  5.1220551  4.869586  4.943840 5.1384776
##         Samples
## Tags       JMS9-P7c  JMS9-P8c
##   497097 -3.2103696 -5.295316
##   20671  -0.3228444 -1.207853
##   27395   4.3396075  4.124644
##   18777   5.7513694  5.142436
##   21399   5.0308985  4.979644
## 16619 more rows ...
## 
## $weights
##           [,1]      [,2]      [,3]      [,4]      [,5]      [,6]      [,7]
## [1,]  1.079413  1.332986 19.826915 20.273253  1.993686  1.395853 20.494977
## [2,]  1.170357  1.456380  4.804866  8.660025  3.612508  2.626870  8.760149
## [3,] 20.219073 25.573792 30.434759 28.528310 31.352260 25.743247 28.722497
## [4,] 26.947557 32.505933 33.583128 33.232125 34.231754 32.354158 33.334340
## [5,] 26.610864 28.501638 33.645479 33.206374 33.573492 31.996623 33.308490
##           [,8]      [,9]
## [1,]  1.107780  1.079413
## [2,]  3.211473  2.541942
## [3,] 21.200072 16.657930
## [4,] 30.348630 24.259801
## [5,] 25.171513 23.573305
## 16619 more rows ...
## 
## $design
##   Basal LP ML laneL006 laneL008
## 1     0  1  0        0        0
## 2     0  0  1        0        0
## 3     1  0  0        0        0
## 4     1  0  0        1        0
## 5     0  0  1        1        0
## 6     0  1  0        1        0
## 7     1  0  0        1        0
## 8     0  0  1        0        1
## 9     0  1  0        0        1
## attr(,"assign")
## [1] 1 1 1 2 2
## attr(,"contrasts")
## attr(,"contrasts")$group
## [1] "contr.treatment"
## 
## attr(,"contrasts")$lane
## [1] "contr.treatment"

vfit <- lmFit(v, design)
vfit <- contrasts.fit(vfit, contrasts=contr.matrix)
efit <- eBayes(vfit)
plotSA(efit, main="Final model: Mean-variance trend")

Si no se tienen más de 3 réplicas se opta por usar un GLM o si se desea usar un modelo específico para RNA-seq.

## Adjust GLM
x <- estimateDisp(x, design)
fit <- glmQLFit(x, design)

## Make contrast matrices
contr.matrix <- makeContrasts(
   BasalvsLP = Basal-LP, 
   BasalvsML = Basal - ML, 
   LPvsML = LP - ML, 
   levels = colnames(design))
contr.matrix

##           Contrasts
## Levels     BasalvsLP BasalvsML LPvsML
##   Basal            1         1      0
##   LP              -1         0      1
##   ML               0        -1     -1
##   laneL006         0         0      0
##   laneL008         0         0      0

## Fit the model
qlf_basal_vs_lp <- glmQLFTest(fit, contrast=contr.matrix[,"BasalvsLP"])
qlf_basal_vs_ml <- glmQLFTest(fit, contrast=contr.matrix[,"BasalvsML"])

## Results
top_basal_vs_lp <- topTags(qlf_basal_vs_lp, n=Inf)
head(top_basal_vs_lp$table)

top_basal_vs_ml <- topTags(qlf_basal_vs_ml, n=Inf)
head(top_basal_vs_ml$table)

## Visualización del ajuste de dispersión
plotQLDisp(fit, main="Dispersion trend under GLM-QL")

## Summary of the results
summary(decideTests(qlf_basal_vs_lp))

##        1*Basal -1*LP
## Down            4613
## NotSig          7032
## Up              4979

summary(decideTests(qlf_basal_vs_ml))

##        1*Basal -1*ML
## Down            4889
## NotSig          6852
## Up              4883

2.7.1 Examining the number of DE genes

La significancia está dada por el threshold al p-value ajustado (5%). Aunque a veces se emplea el valor del log-foldchange como threshold.

## Using adjusted p-value
summary(decideTests(efit))

##        BasalvsLP BasalvsML LPvsML
## Down        4648      4927   3135
## NotSig      7113      7026  10972
## Up          4863      4671   2517

## Using log-foldchange
tfit <- treat(vfit, lfc=1)
dt <- decideTests(tfit)
summary(dt)

##        BasalvsLP BasalvsML LPvsML
## Down        1632      1777    224
## NotSig     12976     12790  16210
## Up          2016      2057    190

## Extract the DEG's shared between both conditions
de.common <- which(dt[,1]!=0 & dt[,2]!=0)
length(de.common)

## [1] 2784

vennDiagram(dt[,1:2], circle.col=c("turquoise", "salmon"))

2.7.1 Examining individual DE genes from top to bottom

Extraer los resultados más significativos.

## CHeck ?topTreat for more details
basal.vs.lp <- topTreat(tfit, coef=1, n=Inf)

basal.vs.ml <- topTreat(tfit, coef=2, n=Inf)
head(basal.vs.lp)

head(basal.vs.ml)

2.7.1 Useful graphical representations of differential expression results

MD plot

plotMD(tfit, column=1, status=dt[,1], main=colnames(tfit)[1], 
       xlim=c(-8,13))

Volcano plot

par(mfrow=c(1,2))

plotVolcano(df = basal.vs.lp,
            pval_col = "P.Value",
            lfc_col = "logFC",
            adjpval_col = "adj.P.Val",
            plot_title = "Basal vs. LP")

## [1] TRUE

plotVolcano(df = basal.vs.ml,
            pval_col = "P.Value",
            lfc_col = "logFC",
            adjpval_col = "adj.P.Val",
            plot_title = "Basal vs. LP")

## [1] TRUE

Heatmap

## Basal vs. LP heatmap
basal.vs.lp.topgenes <- basal.vs.lp$ENTREZID[1:100]
i <- which(v$genes$ENTREZID %in% basal.vs.lp.topgenes)

hmap.palette <- colorRampPalette(c("red", "white", "blue"))
pheatmap::pheatmap(lcpm[i,],
                   scale = "row",
                   cluster_rows = TRUE,
                   labels_row = v$genes$SYMBOL[i],
                   labels_col = group,
                   color = hmap.palette(100),
                   main = "Basal vs. LP")

## Basal vs. ML heatmap
basal.vs.ml.topgenes <- basal.vs.ml$ENTREZID[1:100]
i <- which(v$genes$ENTREZID %in% basal.vs.ml.topgenes)

hmap.palette <- colorRampPalette(c("red", "white", "blue"))
pheatmap::pheatmap(lcpm[i,],
                   scale = "row",
                   cluster_rows = TRUE,
                   labels_row = v$genes$SYMBOL[i],
                   labels_col = group,
                   color = hmap.palette(100),
                   main = "Basal vs. ML")

3. DESeq2

DESeq2 es ubna librería de R que sirve para realizar experimentos de expresión diferencial. Su input es la matriz de conteos crudos resultantes de un experimento de secuenciación de RNA-seq.

• Documentación (DESeq2)
• DOI: https://doi.org/10.1186/s13059-014-0550-8

if (!require("BiocManager", quietly = TRUE))
    install.packages("BiocManager", quiet = TRUE)

if (!require(DESeq2, quietly = TRUE))
    BiocManager::install("DESeq2",quiet = TRUE)

3.1 Load libraries

suppressPackageStartupMessages(library(DESeq2))

3.2 Read files

files <- file.path(tempdir(),"GSE63310/",files)

## Load all the files into a list of dataframes
counts.list <- lapply(files, function(file) {
    data <- read.delim(file, header = TRUE, stringsAsFactors = FALSE)
    data[, c("EntrezID", "Count")]
})
names(counts.list) <- substring(basename(files), 12, nchar(basename(files)) - 4)

## Merge the dataframes into a single counts matrix
counts.mtx <- Reduce(function(x, y) {
    merge(x, y, by = "EntrezID", all = TRUE)
}, counts.list)
rownames(counts.mtx) <- counts.mtx$EntrezID
counts.mtx <- counts.mtx[, -1]
colnames(counts.mtx) <- names(counts.list)

3.2.1 Create metadata

Es neceario que se cree un dataframe con los metadatos del experimento. Aquí una muestra de las condiciones experimentales.

• Group: condiciones experimentales

• Lane: se refiere al batch (lote, ronda de secuenciación)

• lib.size: se trata del total de counts registrados

• norm.factors: factores de normalización

## Create metadata
group <- as.factor(c("LP", "ML", "Basal", "Basal", "ML", "LP",
                     "Basal", "ML", "LP"))
lane <- as.factor(rep(c("L004","L006","L008"), c(3,4,2)))
dds_metadata <- data.frame(group = group,
                           lib.size = colSums(counts.mtx),
                           norm.factors = 1,
                           lane = lane)
dds_metadata

3.3 Create DESeq2 object

Se puede importar de diferentes fuentes. Una de ellas es usando archivos de la abundancia de transcritos con la librería tximport, otra es con una matriz de conteos, otra con archivos de conteos de htseq-count y la última de un objeto de SummarizedExperiment.

## Create DESeq2 object
dds <- DESeqDataSetFromMatrix(countData = counts.mtx,
                              colData = dds_metadata,
                              design = ~ group + lane)
dds

## class: DESeqDataSet 
## dim: 27179 9 
## metadata(1): version
## assays(1): counts
## rownames(27179): 11287 11298 ... 100862405 100862406
## rowData names(0):
## colnames(9): 10_6_5_11 9_6_5_11 ... JMS9-P7c JMS9-P8c
## colData names(4): group lib.size norm.factors lane

3.4 Data pre-filter

Se necesita realizar un pre filtrado de los genes que tuvieron un conteo bajo. Por un lado se reduce el tamaño que ocupa en la memoria el objeto dds e incrementar la velocidad de modelado de los conteos. Como valores standard se ocupa lo siguiente:

• rows con un mínimo de 10 counts.

• Usar el smallest group size el cual es igual al número de muestras.

## We have 9 samples
smallestGroupSize <- 9
keep <- rowSums(counts(dds) >= 10) >= smallestGroupSize
dds <- dds[keep,]
dds

## class: DESeqDataSet 
## dim: 13274 9 
## metadata(1): version
## assays(1): counts
## rownames(13274): 11302 11303 ... 100862307 100862400
## rowData names(0):
## colnames(9): 10_6_5_11 9_6_5_11 ... JMS9-P7c JMS9-P8c
## colData names(4): group lib.size norm.factors lane

3.4.1 Refactor metadata

## Put the reference in the first position of the levels
dds$group <- factor(dds$group, levels = c("Basal","LP","ML"))

## Or just relevel them specifying the reference
dds$group <- relevel(dds$group, ref = "Basal")

Adicionalmente se pueden colapsar todas las réplicas en columnas individuales con la función DESeq2::collapseReplicates() ESTO SOLO PARA RÉPLICAS TÉCNICAS, NO PARA BIOLÓGICAS (una réplica técnica implica múltiples corridas de secuenciación de la misma librería).

3.4.2 Add gene names

genes <- AnnotationDbi::select(Mus.musculus,
                        keys = rownames(dds),
                        columns = c("SYMBOL","TXCHROM"),
                        keytype = "ENTREZID")

## Remove duplicates
genes <- genes[!duplicated(genes$ENTREZID),]
rownames(genes) <- genes$ENTREZID

## Add genes
rowData(dds) <- genes

3.5 Differential Expression Analysis

Así de simple.

dds <- DESeq(dds)
res <- results(dds)
res

## log2 fold change (MLE): lane L008 vs L004 
## Wald test p-value: lane L008 vs L004 
## DataFrame with 13274 rows and 6 columns
##            baseMean log2FoldChange     lfcSE      stat    pvalue      padj
##           <numeric>      <numeric> <numeric> <numeric> <numeric> <numeric>
## 11302      524.8568      0.2025970  0.427536  0.473871  0.635592  0.914339
## 11303     5688.7471      0.1132199  0.312441  0.362372  0.717074  0.935875
## 11304       43.6222      0.1260263  0.510059  0.247082  0.804845  0.960693
## 11305     3258.0083     -0.0949987  0.168277 -0.564539  0.572387  0.893517
## 11306     1850.6776      0.0372121  0.181120  0.205455  0.837216  0.967711
## ...             ...            ...       ...       ...       ...       ...
## 100862251 1084.2679       0.157915  0.199627  0.791051  0.428914  0.842320
## 100862257 3711.2564      -0.180482  0.157793 -1.143787  0.252712  0.729788
## 100862258   67.0328      -0.518082  0.381632 -1.357545  0.174608  0.646332
## 100862307   51.7442       0.321677  0.394997  0.814379  0.415428  0.834233
## 100862400   88.7577       0.359888  0.497736  0.723049  0.469650  0.857534

Si quisiéramos extraer los genes diferencialmente expresados filtrando por el adjusted p-value entonces se le debe de agregar el parámetro alpha = <float> a la función results(). Por ejemplo, si queremos los DEG’s que son menores o iguales a 0.05 entonces se ve así:

res05 <- results(dds, alpha = 0.05)
summary(res05)

## 
## out of 13274 with nonzero total read count
## adjusted p-value < 0.05
## LFC > 0 (up)       : 112, 0.84%
## LFC < 0 (down)     : 196, 1.5%
## outliers [1]       : 0, 0%
## low counts [2]     : 0, 0%
## (mean count < 17)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results

sum(res05$padj < 0.05, na.rm=TRUE)

## [1] 308

Además, se pueden revisar las comparaciones hechas:

resultsNames(dds)

## [1] "Intercept"         "group_LP_vs_Basal" "group_ML_vs_Basal"
## [4] "lane_L006_vs_L004" "lane_L008_vs_L004"

3.5.1 Exploring results

Para la visualización se pueden “encoger” los resultados de acuerdo al efecto del tamaño. El estimador por default es apeglm el cual estima con una tendencia al 0, lo que promueve la estabilidad de los estimados en los cuales hay una variabilidad alta o tamaños de muestra pequeños.

resLFC.LPBas <- lfcShrink(dds,coef = 2,type = "apeglm")
resLFC.MLBas <- lfcShrink(dds,coef = 3,type = "apeglm")

MA-plot

plotMA(res05, ylim=c(-2,2))

plotMA(resLFC.LPBas, ylim=c(-2,2))

plotMA(resLFC.MLBas, ylim=c(-2,2))

Get dispersion factors

Aquí se obtienen los factores de dispersión que es similar a lo que hace la función voom en EdgeR. Aquí se usa la función para Variance stabilizing transformations (vst) la cual transforma a los datos en una escala log2 la cual ha sido normalizada con respecto al tamaño de la librería u otros factores de normalización. La finalidad de esto es eliminar la dependencia de la varianza en la media (si se trabaja con niveles de expresión, entonces genes con baja expresión tendrán una desviación estándar mayor que los genes de alta expresión debido al tipo de distribución de los datos de RNA-seq). El parámetro blind es TRUE por default porque este re-estima las dispersiones empleando solo un intercepto. Esto no es adecuado si se espera que la mayoría de los genes tengan grandes diferencias en sus conteos (lo que puede causar grandes estimados de dispersión por el diseño experimental y no por una diferencia biológica).

## Variance stabilizing transformation
vsd <- vst(dds,blind = FALSE)

## Plot
vsn::meanSdPlot(assay(vsd))

3.5.2 Extract the top DEG’s from top to bottom

Una vez que ya se tiene la normalización con vsd, se puede proceder a extraer los DEG’s de interés para cada una de las comparaciones.

## Extract the results for each condition
res_LP_vs_Basal <- results(dds, contrast = c("group", "LP", "Basal"))
res_ML_vs_Basal <- results(dds, contrast = c("group", "ML", "Basal"))

## Add gene names
res_LP_vs_Basal$SYMBOL <- rowData(dds[rownames(res_LP_vs_Basal)])$SYMBOL
res_ML_vs_Basal$SYMBOL <- rowData(dds[rownames(res_ML_vs_Basal)])$SYMBOL
res_LP_vs_Basal$ENTREZID <- rowData(dds[rownames(res_LP_vs_Basal)])$ENTREZID
res_ML_vs_Basal$ENTREZID <- rowData(dds[rownames(res_ML_vs_Basal)])$ENTREZID

## Filter the results
res_LP_vs_Basal_sig <- res_LP_vs_Basal[which(res_LP_vs_Basal$padj < 0.05), ]
res_LP_vs_Basal_sig <- res_LP_vs_Basal_sig[order(res_LP_vs_Basal_sig$padj, decreasing = FALSE), ]

res_ML_vs_Basal_sig <- res_ML_vs_Basal[which(res_ML_vs_Basal$padj < 0.05), ]
res_ML_vs_Basal_sig <- res_ML_vs_Basal_sig[order(res_ML_vs_Basal_sig$padj, decreasing = FALSE), ]

Volcano plot

par(mfrow=c(1,2))

plotVolcano(df = res_LP_vs_Basal_sig,
            pval_col = "pvalue",
            lfc_col = "log2FoldChange",
            adjpval_col = "padj",
            plot_title = "Basal vs. LP")

## [1] TRUE

plotVolcano(df = res_ML_vs_Basal_sig,
            pval_col = "pvalue",
            lfc_col = "log2FoldChange",
            adjpval_col = "padj",
            plot_title = "Basal vs. ML")

## [1] TRUE

Heatmap

## Basal vs. LP heatmap
basal.vs.lp.topgenes.deseq <- res_LP_vs_Basal_sig$ENTREZID[1:100]
i <- which(rowData(dds)$ENTREZID %in% basal.vs.lp.topgenes.deseq)

hmap.palette <- colorRampPalette(c("red", "white", "blue"))
pheatmap::pheatmap(assay(vsd)[i,],
                   scale = "row",
                   cluster_rows = TRUE,
                   labels_row = rowData(dds[i])$SYMBOL,
                   labels_col = group,
                   color = hmap.palette(100),
                   main = "Basal vs. LP")

## Basal vs. ML heatmap
basal.vs.ml.topgenes.deseq <- res_ML_vs_Basal_sig$ENTREZID[1:100]
i <- which(rowData(dds)$ENTREZID %in% basal.vs.ml.topgenes.deseq)

hmap.palette <- colorRampPalette(c("red", "white", "blue"))
pheatmap::pheatmap(assay(vsd)[i,],
                   scale = "row",
                   cluster_rows = TRUE,
                   labels_row = rowData(dds[i])$SYMBOL,
                   labels_col = group,
                   color = hmap.palette(100),
                   main = "Basal vs. ML")

4. pyDESeq2

Este es una librería de Pyhton la cual pretende implementar la librería original de R DESeq2 para el análisis de expresión diferencial con datos de bulk RNA-seq. Como es una implementación desde 0, hay algunos cambios importantes en términos de las funciones y de los valores que estas regresan.

• Documentación (pyDESeq2)
• DOI: https://doi.org/10.1093/bioinformatics/btad547

pip install pydeseq2

## Using conda
conda env create -n pydeseq2
conda activate pydeseq2
pip install pydeseq2

4.1 Load libraries

import re
import numpy as np
import seaborn as sns
import pandas as pd
import matplotlib.pyplot as plt

# %pip install sanbomics
from sanbomics.plots import volcano
from sanbomics.tools import id_map

# %pip install pydeseq2
from pydeseq2.dds import DeseqDataSet
from pydeseq2.default_inference import DefaultInference
from pydeseq2.ds import DeseqStats

# %pip install scanpy
import scanpy as sc

4.2 Read files

## Import R files variable to python
files = r.files

## Empty df
dfs = []

## Read files
for file in files:
    df = pd.read_csv(file, sep="\t",usecols = ["EntrezID","Count"])
    column_name = file.split("/")[-1].replace(".txt","")
    column_name = re.sub(r"^GSM15455\d{2}_", "", column_name)
    df.rename(columns={"Count": column_name}, inplace=True)
    
    dfs.append(df)

## Create the matrix
merged_matrix = pd.concat(
    [df.set_index("EntrezID") for df in dfs], axis=1
).fillna(0)

4.2.1 Create metadata

metadata = r.dds_metadata

4.2.2 Add genes

genes = r.genes

4.2.3 Pre-filter counts

## We have 9 samples
smallest_group_size = 9
keep = np.sum(merged_matrix >= 10, axis=1) >= smallest_group_size
merged_matrix = merged_matrix[keep].copy()

4.3 Create pyDESeq2 object

## Create DESeq2 object
inference = DefaultInference(n_cpus = 8)
dds = DeseqDataSet(
    counts = merged_matrix.T,
    metadata = metadata,
    design_factors = ['group', 'lane'],
    inference = inference
)
dds.var = genes

4.4 Differential Expression Analysis

## Run DESeq2
dds.deseq2()

## Fitting size factors...
## ... done in 0.00 seconds.
## 
## Fitting dispersions...
## ... done in 1.28 seconds.
## 
## Fitting dispersion trend curve...
## ... done in 0.15 seconds.
## 
## Fitting MAP dispersions...
## ... done in 1.35 seconds.
## 
## Fitting LFCs...
## ... done in 0.69 seconds.
## 
## Calculating cook's distance...
## ... done in 0.00 seconds.
## 
## Replacing 0 outlier genes.

## Retrieve stats
result = DeseqStats(dds,inference = inference)
basal_vs_lp = DeseqStats(dds, contrast=["group", "Basal", "LP"], inference=inference)
basal_vs_ml = DeseqStats(dds, contrast=["group", "Basal", "ML"], inference=inference)

## Extract the results
result.summary()

## Running Wald tests...
## ... done in 0.38 seconds.
## 
## Log2 fold change & Wald test p-value: lane L006 vs L004
##               baseMean  log2FoldChange     lfcSE      stat    pvalue      padj
## 11302       782.020723       -0.123189  0.202669 -0.607834  0.543297  0.875635
## 11303      1651.935164        0.039150  0.159773  0.245032  0.806431  0.959511
## 11304      1344.515458        0.006769  0.135751  0.049865  0.960230  0.992764
## 11305       295.742823       -0.399928  0.211048 -1.894967  0.058097  0.416339
## 11306      2356.502410        0.055302  0.137805  0.401304  0.688196  0.926800
## ...                ...             ...       ...       ...       ...       ...
## 100862251   838.908288        0.297839  0.232688  1.279994  0.200547  0.660890
## 100862257  1041.025804       -0.139267  0.196478 -0.708816  0.478438  0.843321
## 100862258  1534.037447       -0.227660  0.191033 -1.191730  0.233367  0.692072
## 100862307   369.536157       -0.355791  0.216680 -1.642006  0.100589  0.523408
## 100862400   855.267041       -0.344427  0.208666 -1.650609  0.098818  0.520729
## 
## [13274 rows x 6 columns]

res = result.results_df
basal_vs_lp.summary()

## Running Wald tests...
## ... done in 0.34 seconds.
## 
## Log2 fold change & Wald test p-value: group Basal vs LP
##               baseMean  log2FoldChange  ...        pvalue          padj
## 11302       782.020723        0.523221  ...  2.276141e-02  4.015082e-02
## 11303      1651.935164        0.319091  ...  7.795831e-02  1.195217e-01
## 11304      1344.515458        0.397116  ...  9.764703e-03  1.865793e-02
## 11305       295.742823       -5.521131  ...  1.947981e-77  1.213967e-75
## 11306      2356.502410       -0.791075  ...  3.972588e-07  1.463156e-06
## ...                ...             ...  ...           ...           ...
## 100862251   838.908288        0.128801  ...  6.241817e-01  6.960756e-01
## 100862257  1041.025804        0.187428  ...  4.008521e-01  4.868171e-01
## 100862258  1534.037447        0.143802  ...  5.059068e-01  5.888641e-01
## 100862307   369.536157        0.047139  ...  8.487170e-01  8.857512e-01
## 100862400   855.267041       -0.072764  ...  7.580740e-01  8.120299e-01
## 
## [13274 rows x 6 columns]

basal_vs_lp_res = basal_vs_lp.results_df
basal_vs_ml.summary()

## Running Wald tests...
## ... done in 0.35 seconds.
## 
## Log2 fold change & Wald test p-value: group Basal vs ML
##               baseMean  log2FoldChange  ...        pvalue          padj
## 11302       782.020723        0.125763  ...  5.822965e-01  6.564249e-01
## 11303      1651.935164       -0.310793  ...  8.461658e-02  1.275209e-01
## 11304      1344.515458        0.460100  ...  2.679736e-03  5.633643e-03
## 11305       295.742823       -5.023438  ...  2.219609e-64  8.299461e-63
## 11306      2356.502410       -0.979360  ...  3.187953e-10  1.490556e-09
## ...                ...             ...  ...           ...           ...
## 100862251   838.908288        0.708760  ...  7.059043e-03  1.372517e-02
## 100862257  1041.025804       -1.106419  ...  5.835175e-07  2.015512e-06
## 100862258  1534.037447        0.310210  ...  1.509529e-01  2.109996e-01
## 100862307   369.536157       -0.459937  ...  6.060102e-02  9.521993e-02
## 100862400   855.267041        0.166803  ...  4.797587e-01  5.609369e-01
## 
## [13274 rows x 6 columns]

basal_vs_ml_res = basal_vs_ml.results_df

4.4.1 Add gene names

## Add the genes information
res = res.merge(dds.var,left_index=True,right_index=True)
basal_vs_lp_res = basal_vs_lp_res.merge(dds.var,left_index=True,right_index=True)
basal_vs_ml_res = basal_vs_ml_res.merge(dds.var,left_index=True,right_index=True)

## Sort by adjusted p-value
basal_vs_lp_res = basal_vs_lp_res.sort_values(by='padj', ascending=True)
basal_vs_ml_res = basal_vs_ml_res.sort_values(by='padj', ascending=True)

4.4.2 Exploring results

MA-plot

## MA-plot
plt.figure(figsize=(8, 6))
sns.scatterplot(
    x=res["baseMean"],
    y=res["log2FoldChange"],
    hue=res["padj"] < 0.05,  # Highlight significant points
    palette={True: "blue", False: "gray"},
    alpha=0.6,
)
plt.axhline(0, color="black", linestyle="--", linewidth=1)
plt.xscale("log")
plt.xlabel("Mean Normalized Counts")
plt.ylabel("Log2 Fold Change")
plt.title("MA Plot")
plt.legend(title="Significant", labels=["Yes", "No"])
plt.tight_layout()
plt.show()

PCA

sc.tl.pca(dds)

## Group PCA
sc.pl.pca(dds, color='group', size=200, annotate_var_explained=True,legend_loc='best', show=True)

## Lane PCA
sc.pl.pca(dds, color='lane', size=200, annotate_var_explained=True,legend_loc='best', show=True)

Volcano plot

## Basal vs. LP volcano plot
volcano(
    basal_vs_lp_res,
    symbol = 'SYMBOL',
    to_label= 30,
    pval_thresh = 0.05,
    log2fc_thresh = 2.5)

## 0s encountered for p value, imputing 1e-323
## impute your own value if you want to avoid this

    
## Basal vs. LP volcano plot
volcano(
    basal_vs_lp_res,
    symbol = 'SYMBOL',
    to_label= 30,
    pval_thresh = 0.05,
    log2fc_thresh = 2.5)

## 0s encountered for p value, imputing 1e-323
## impute your own value if you want to avoid this

Heatmap

## Extract the top genes
basal_vs_lp_topgenes = basal_vs_lp_res.ENTREZID[0:100]
basal_vs_ml_topgenes = basal_vs_ml_res.ENTREZID[0:100]

## Add log transformation counts
dds.layers['log1p'] = np.log1p(dds.layers['normed_counts'])

## Get index from significant genes
basalvslp_sigs = dds[:, basal_vs_lp_topgenes]
basalvsml_sigs = dds[:, basal_vs_ml_topgenes]

## Create counts matrices
grapher_blp = pd.DataFrame(
    basalvslp_sigs.layers['log1p'].T,
    index=basalvslp_sigs.var_names,
    columns=basalvslp_sigs.obs_names)

grapher_bml = pd.DataFrame(
    basalvsml_sigs.layers['log1p'].T,
    index=basalvsml_sigs.var_names,
    columns=basalvsml_sigs.obs_names)

## Define color palette
hmap_palette = sns.diverging_palette(10, 250, as_cmap=True)

## Heatmap for Basal vs. LP
sns.clustermap(grapher_blp,
               z_score=0,
               cmap=hmap_palette,
               col_cluster=True,
               row_cluster=True,
               xticklabels=metadata.group[grapher_bml.columns].to_list(),
               yticklabels=basal_vs_lp_res.SYMBOL[basal_vs_lp_topgenes].to_list()).fig.suptitle('Basal vs. LP') 

plt.show()

## Heatmap for Basal vs. ML
sns.clustermap(grapher_bml,
               z_score=0,
               cmap=hmap_palette,
               col_cluster=True,
               row_cluster=True,
               xticklabels=metadata.group[grapher_bml.columns],
               yticklabels=basal_vs_ml_res.SYMBOL[basal_vs_ml_topgenes].to_list()).fig.suptitle('Basal vs. ML') 

plt.show()

5. Comparación

5.1 edgeR

• Método principal:
  Utiliza modelos lineales generalizados (GLM) basados en la distribución binomial negativa.

• Normalización:
  Realiza la normalización utilizando el método de la media recortada ponderada de la razón (TMM). Este método ajusta los efectos de composición entre muestras.

• Estadísticas diferenciales:
  Estima la dispersión biológica para cada gen y utiliza un enfoque Bayesiano empírico para estabilizar estas estimaciones, mejorando la precisión en experimentos con pocos replicados.

• Casos de uso:

  • Ideal para datos con pocas réplicas biológicas.

  • Muy robusto frente a variaciones en la profundidad de secuenciación.

5.2 DESeq2

• Método principal:
  Basado en la distribución binomial negativa, modela los datos de conteo utilizando un enfoque de regresión GLM, similar a edgeR. Sin embargo, DESeq2 implementa un procedimiento más detallado para la estimación de dispersión.

• Normalización:
  Utiliza el método de razón geométrica para calcular factores de escala, que corrigen las diferencias de tamaño de biblioteca entre muestras.

• Estadísticas diferenciales:

  • Estima la dispersión utilizando un modelo empírico y luego ajusta estos valores utilizando un enfoque Bayesiano para estabilizar las varianzas.

  • Implementa un procedimiento de “shrinkage” para los coeficientes del modelo, lo que mejora la robustez en genes de baja expresión.

• Casos de uso:

  • Preferido para datos con un número moderado a grande de réplicas.

  • Mejora la detección de genes diferencialmente expresados con bajos conteos.

5.3 pyDESeq2

• Método principal:
  Es una implementación en Python basada en DESeq2. Internamente utiliza las mismas bases teóricas y estadísticas, ya que está diseñado para ofrecer la misma funcionalidad en un entorno Python.

• Normalización y estadística:

  • Reproduce la normalización basada en razones geométricas y los modelos de regresión GLM con dispersión ajustada Bayesiana.

  • La diferencia clave está en que facilita la integración con bibliotecas y pipelines en Python, como pandas y scikit-learn.

• Casos de uso:

  • Ideal para usuarios que trabajan en ecosistemas Python y no desean depender de R.

  • Útil en pipelines integrados con otras herramientas de aprendizaje automático o análisis avanzado.

5.4 Resumen

Aspecto	edgeR	DESeq2	pyDESeq2
Base estadística	Binomial negativa + GLM	Binomial negativa + GLM	Igual que DESeq2
Normalización	TMM	Razón geométrica	Igual que DESeq2
Ajuste Bayesiano	Sí (dispersión)	Sí (dispersión y shrinkage)	Igual que DESeq2
Facilidad de uso	R	R	Python
Casos recomendados	Pocas réplicas	Réplicas moderadas o altas	Usuarios de Python

5.5 Resultados obtenidos

Hay que recordar que pyDESeq2 es una librería emergente que trata de implementar la librería original de R DESeq2 en Python. Los resultados muestran una similitud de 86 y 84 genes cuando se comparan las librerías de R DESeq2 y edgeR mientras que se encuentra sólo 1 gen compartido con los resultados de pyDESeq2. Esto mayormente puede deberse a la implementación dada. Sin embargo, no hay que descartar que esten ahí hasta una posterior confirmación (probablemente sólo estén rankeados diferente). Por eso la recomendación es que se sigan usando las librerías de R hasta que se logre hacer una transición completa a Python en términos del análisis transcriptómico.

Reproducibilidad

R

## ─ Session info ───────────────────────────────────────────────────────────────
##  setting  value
##  version  R version 4.4.2 (2024-10-31 ucrt)
##  os       Windows 11 x64 (build 26100)
##  system   x86_64, mingw32
##  ui       RTerm
##  language (EN)
##  collate  English_United States.utf8
##  ctype    English_United States.utf8
##  tz       America/Mexico_City
##  date     2025-01-21
##  pandoc   3.2 @ C:/Program Files/RStudio/resources/app/bin/quarto/bin/tools/ (via rmarkdown)
## 
## ─ Packages ───────────────────────────────────────────────────────────────────
##  package                            * version    date (UTC) lib source
##  abind                                1.4-8      2024-09-12 [1] CRAN (R 4.4.1)
##  affy                                 1.84.0     2024-10-29 [1] Bioconductor 3.20 (R 4.4.1)
##  affyio                               1.76.0     2024-10-29 [1] Bioconductor 3.20 (R 4.4.1)
##  AnnotationDbi                      * 1.68.0     2024-10-29 [1] Bioconductor 3.20 (R 4.4.1)
##  apeglm                               1.28.0     2024-10-29 [1] Bioconductor 3.20 (R 4.4.1)
##  bbmle                                1.0.25.1   2023-12-09 [1] CRAN (R 4.4.2)
##  bdsmatrix                            1.3-7      2024-03-02 [1] CRAN (R 4.4.0)
##  Biobase                            * 2.66.0     2024-10-29 [1] Bioconductor 3.20 (R 4.4.1)
##  BiocFileCache                        2.14.0     2024-10-29 [1] Bioconductor 3.20 (R 4.4.1)
##  BiocGenerics                       * 0.52.0     2024-10-29 [1] Bioconductor 3.20 (R 4.4.1)
##  BiocIO                               1.16.0     2024-10-29 [1] Bioconductor 3.20 (R 4.4.1)
##  BiocManager                          1.30.25    2024-08-28 [1] CRAN (R 4.4.2)
##  BiocParallel                         1.40.0     2024-10-29 [1] Bioconductor 3.20 (R 4.4.1)
##  biomaRt                              2.62.0     2024-10-29 [1] Bioconductor 3.20 (R 4.4.1)
##  Biostrings                           2.74.1     2024-12-16 [1] Bioconductor 3.20 (R 4.4.2)
##  bit                                  4.5.0.1    2024-12-03 [1] CRAN (R 4.4.2)
##  bit64                                4.5.2      2024-09-22 [1] CRAN (R 4.4.2)
##  bitops                               1.0-9      2024-10-03 [1] CRAN (R 4.4.1)
##  blob                                 1.2.4      2023-03-17 [1] CRAN (R 4.4.2)
##  cachem                               1.1.0      2024-05-16 [1] CRAN (R 4.4.2)
##  cli                                  3.6.3      2024-06-21 [1] CRAN (R 4.4.2)
##  coda                                 0.19-4.1   2024-01-31 [1] CRAN (R 4.4.2)
##  codetools                            0.2-20     2024-03-31 [2] CRAN (R 4.4.2)
##  colorspace                           2.1-1      2024-07-26 [1] CRAN (R 4.4.2)
##  crayon                               1.5.3      2024-06-20 [1] CRAN (R 4.4.2)
##  curl                                 6.0.1      2024-11-14 [1] CRAN (R 4.4.2)
##  DBI                                  1.2.3      2024-06-02 [1] CRAN (R 4.4.2)
##  dbplyr                               2.5.0      2024-03-19 [1] CRAN (R 4.4.2)
##  DelayedArray                         0.32.0     2024-10-29 [1] Bioconductor 3.20 (R 4.4.1)
##  DESeq2                             * 1.46.0     2024-10-29 [1] Bioconductor 3.20 (R 4.4.1)
##  digest                               0.6.37     2024-08-19 [1] CRAN (R 4.4.2)
##  dplyr                                1.1.4      2023-11-17 [1] CRAN (R 4.4.2)
##  edgeR                              * 4.4.1      2024-12-02 [1] Bioconductor 3.20 (R 4.4.2)
##  emdbook                              1.3.13     2023-07-03 [1] CRAN (R 4.4.2)
##  evaluate                             1.0.1      2024-10-10 [1] CRAN (R 4.4.2)
##  farver                               2.1.2      2024-05-13 [1] CRAN (R 4.4.2)
##  fastmap                              1.2.0      2024-05-15 [1] CRAN (R 4.4.2)
##  filelock                             1.0.3      2023-12-11 [1] CRAN (R 4.4.2)
##  generics                             0.1.3      2022-07-05 [1] CRAN (R 4.4.2)
##  GenomeInfoDb                       * 1.42.1     2024-11-28 [1] Bioconductor 3.20 (R 4.4.2)
##  GenomeInfoDbData                     1.2.13     2024-12-22 [1] Bioconductor
##  GenomicAlignments                    1.42.0     2024-10-29 [1] Bioconductor 3.20 (R 4.4.1)
##  GenomicFeatures                    * 1.58.0     2024-10-29 [1] Bioconductor 3.20 (R 4.4.1)
##  GenomicRanges                      * 1.58.0     2024-10-29 [1] Bioconductor 3.20 (R 4.4.1)
##  ggplot2                              3.5.1      2024-04-23 [1] CRAN (R 4.4.2)
##  ggVennDiagram                      * 1.5.2      2024-02-20 [1] CRAN (R 4.4.2)
##  Glimma                             * 2.16.0     2024-10-29 [1] Bioconductor 3.20 (R 4.4.1)
##  glue                                 1.8.0      2024-09-30 [1] CRAN (R 4.4.2)
##  GO.db                              * 3.20.0     2024-12-22 [1] Bioconductor
##  graph                                1.84.0     2024-10-29 [1] Bioconductor 3.20 (R 4.4.1)
##  gtable                               0.3.6      2024-10-25 [1] CRAN (R 4.4.2)
##  hexbin                               1.28.5     2024-11-13 [1] CRAN (R 4.4.2)
##  hms                                  1.1.3      2023-03-21 [1] CRAN (R 4.4.2)
##  htmltools                            0.5.8.1    2024-04-04 [1] CRAN (R 4.4.2)
##  htmlwidgets                          1.6.4      2023-12-06 [1] CRAN (R 4.4.2)
##  httr                                 1.4.7      2023-08-15 [1] CRAN (R 4.4.2)
##  httr2                                1.0.7      2024-11-26 [1] CRAN (R 4.4.2)
##  IRanges                            * 2.40.1     2024-12-05 [1] Bioconductor 3.20 (R 4.4.2)
##  jsonlite                             1.8.9      2024-09-20 [1] CRAN (R 4.4.2)
##  KEGGREST                             1.46.0     2024-10-29 [1] Bioconductor 3.20 (R 4.4.1)
##  KernSmooth                           2.23-24    2024-05-17 [2] CRAN (R 4.4.2)
##  knitr                                1.49       2024-11-08 [1] CRAN (R 4.4.2)
##  labeling                             0.4.3      2023-08-29 [1] CRAN (R 4.4.0)
##  lattice                              0.22-6     2024-03-20 [2] CRAN (R 4.4.2)
##  lifecycle                            1.0.4      2023-11-07 [1] CRAN (R 4.4.2)
##  limma                              * 3.62.1     2024-11-03 [1] Bioconductor 3.20 (R 4.4.1)
##  locfit                               1.5-9.10   2024-06-24 [1] CRAN (R 4.4.2)
##  magrittr                             2.0.3      2022-03-30 [1] CRAN (R 4.4.2)
##  MASS                                 7.3-61     2024-06-13 [2] CRAN (R 4.4.2)
##  Matrix                               1.7-1      2024-10-18 [2] CRAN (R 4.4.2)
##  MatrixGenerics                     * 1.18.0     2024-10-29 [1] Bioconductor 3.20 (R 4.4.1)
##  matrixStats                        * 1.4.1      2024-09-08 [1] CRAN (R 4.4.2)
##  memoise                              2.0.1      2021-11-26 [1] CRAN (R 4.4.2)
##  munsell                              0.5.1      2024-04-01 [1] CRAN (R 4.4.2)
##  Mus.musculus                       * 1.3.1      2024-12-27 [1] Bioconductor
##  mvtnorm                              1.3-2      2024-11-04 [1] CRAN (R 4.4.2)
##  numDeriv                             2016.8-1.1 2019-06-06 [1] CRAN (R 4.4.0)
##  org.Mm.eg.db                       * 3.20.0     2024-12-22 [1] Bioconductor
##  OrganismDbi                        * 1.48.0     2024-10-29 [1] Bioconductor 3.20 (R 4.4.1)
##  pheatmap                             1.0.12     2019-01-04 [1] CRAN (R 4.4.2)
##  pillar                               1.10.0     2024-12-17 [1] CRAN (R 4.4.2)
##  pkgconfig                            2.0.3      2019-09-22 [1] CRAN (R 4.4.2)
##  plyr                                 1.8.9      2023-10-02 [1] CRAN (R 4.4.2)
##  png                                  0.1-8      2022-11-29 [1] CRAN (R 4.4.0)
##  preprocessCore                       1.68.0     2024-10-29 [1] Bioconductor 3.20 (R 4.4.1)
##  prettyunits                          1.2.0      2023-09-24 [1] CRAN (R 4.4.2)
##  progress                             1.2.3      2023-12-06 [1] CRAN (R 4.4.2)
##  R.methodsS3                          1.8.2      2022-06-13 [1] CRAN (R 4.4.0)
##  R.oo                                 1.27.0     2024-11-01 [1] CRAN (R 4.4.1)
##  R.utils                              2.12.3     2023-11-18 [1] CRAN (R 4.4.2)
##  R6                                   2.5.1      2021-08-19 [1] CRAN (R 4.4.2)
##  rappdirs                             0.3.3      2021-01-31 [1] CRAN (R 4.4.2)
##  RBGL                                 1.82.0     2024-10-29 [1] Bioconductor 3.20 (R 4.4.1)
##  RColorBrewer                         1.1-3      2022-04-03 [1] CRAN (R 4.4.0)
##  Rcpp                                 1.0.13-1   2024-11-02 [1] CRAN (R 4.4.2)
##  RCurl                                1.98-1.16  2024-07-11 [1] CRAN (R 4.4.1)
##  restfulr                             0.0.15     2022-06-16 [1] CRAN (R 4.4.2)
##  reticulate                         * 1.40.0     2024-11-15 [1] CRAN (R 4.4.2)
##  rjson                                0.2.23     2024-09-16 [1] CRAN (R 4.4.1)
##  rlang                                1.1.4      2024-06-04 [1] CRAN (R 4.4.2)
##  rmarkdown                            2.29       2024-11-04 [1] CRAN (R 4.4.2)
##  Rsamtools                            2.22.0     2024-10-29 [1] Bioconductor 3.20 (R 4.4.1)
##  RSQLite                              2.3.9      2024-12-03 [1] CRAN (R 4.4.2)
##  rstudioapi                           0.17.1     2024-10-22 [1] CRAN (R 4.4.2)
##  rtracklayer                          1.66.0     2024-10-29 [1] Bioconductor 3.20 (R 4.4.1)
##  S4Arrays                             1.6.0      2024-10-29 [1] Bioconductor 3.20 (R 4.4.1)
##  S4Vectors                          * 0.44.0     2024-10-29 [1] Bioconductor 3.20 (R 4.4.1)
##  scales                               1.3.0      2023-11-28 [1] CRAN (R 4.4.2)
##  sessioninfo                          1.2.2      2021-12-06 [1] CRAN (R 4.4.2)
##  SparseArray                          1.6.0      2024-10-29 [1] Bioconductor 3.20 (R 4.4.1)
##  statmod                              1.5.0      2023-01-06 [1] CRAN (R 4.4.2)
##  stringi                              1.8.4      2024-05-06 [1] CRAN (R 4.4.0)
##  stringr                              1.5.1      2023-11-14 [1] CRAN (R 4.4.2)
##  SummarizedExperiment               * 1.36.0     2024-10-29 [1] Bioconductor 3.20 (R 4.4.1)
##  tibble                               3.2.1      2023-03-20 [1] CRAN (R 4.4.2)
##  tidyselect                           1.2.1      2024-03-11 [1] CRAN (R 4.4.2)
##  TxDb.Mmusculus.UCSC.mm10.knownGene * 3.10.0     2024-12-22 [1] Bioconductor
##  txdbmaker                            1.2.1      2024-11-25 [1] Bioconductor 3.20 (R 4.4.2)
##  UCSC.utils                           1.2.0      2024-10-29 [1] Bioconductor 3.20 (R 4.4.1)
##  vctrs                                0.6.5      2023-12-01 [1] CRAN (R 4.4.2)
##  vsn                                  3.74.0     2024-10-29 [1] Bioconductor 3.20 (R 4.4.1)
##  withr                                3.0.2      2024-10-28 [1] CRAN (R 4.4.2)
##  xfun                                 0.49       2024-10-31 [1] CRAN (R 4.4.2)
##  XML                                  3.99-0.17  2024-06-25 [1] CRAN (R 4.4.1)
##  xml2                                 1.3.6      2023-12-04 [1] CRAN (R 4.4.2)
##  XVector                              0.46.0     2024-10-29 [1] Bioconductor 3.20 (R 4.4.1)
##  yaml                                 2.3.10     2024-07-26 [1] CRAN (R 4.4.1)
##  zlibbioc                             1.52.0     2024-10-29 [1] Bioconductor 3.20 (R 4.4.1)
## 
##  [1] C:/Users/chomb/AppData/Local/R/win-library/4.4
##  [2] C:/Program Files/R/R-4.4.2/library
## 
## ─ Python configuration ───────────────────────────────────────────────────────
##  python:         C:/Users/chomb/AppData/Local/Programs/Python/Python312/python.exe
##  libpython:      C:/Users/chomb/AppData/Local/Programs/Python/Python312/python312.dll
##  pythonhome:     C:/Users/chomb/AppData/Local/Programs/Python/Python312
##  version:        3.12.4 (tags/v3.12.4:8e8a4ba, Jun  6 2024, 19:30:16) [MSC v.1940 64 bit (AMD64)]
##  Architecture:   64bit
##  numpy:          C:/Users/chomb/AppData/Local/Programs/Python/Python312/Lib/site-packages/numpy
##  numpy_version:  2.0.2
##  
##  NOTE: Python version was forced by use_python() function
## 
## ──────────────────────────────────────────────────────────────────────────────

Python

## -----
## anndata             0.11.1
## matplotlib          3.10.0
## numpy               2.0.2
## pandas              2.2.3
## pydeseq2            0.4.12
## sanbomics           NA
## scanpy              1.10.4
## seaborn             0.13.2
## session_info        1.0.0
## -----
## Python 3.12.4 (tags/v3.12.4:8e8a4ba, Jun  6 2024, 19:30:16) [MSC v.1940 64 bit (AMD64)]
## Windows-11-10.0.26100-SP0
## -----
## Session information updated at 2025-01-21 14:50